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in the Hubbard Model on the Stacked Square Lattice: 
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Using variational cluster approach, we study influence of frustration and dimensionality on mag- 
netic properties in the ground state of Hubbard model on a stacked square lattice in the large U 
region (U/t=10), by changing the next-nearest-neighbor hopping, t', and the interlayer hopping, t±. 
For small t± < t± with i^/i~0.44, antiferromagnetic long-range order appears at small t' < t' ci , 
and collinear magnetic long-range order at large t' > t' C2 . They are separated by a paramagnetic 
Mott insulating state which appears in the parameter region £' < t' < t' C2 . For large t± > t* ± , the 
paramagnetic Mott insulating state disappears and a direct transition between antiferromagnetic 
state and collinear magnetic state occurs at t' — t' C3 . We also find that the transition from the 
antiferromagnetic state to the paramagnetic Mott insulating state is of the second-order and that 
from the collinear magnetic state to the paramagnetic Mott insulating or the antiferromagnetic state 
is of the first-order. 

PACS numbers: 71.10.Fd, 75.10.Jm, 75.40.Mg 
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I. INTRODUCTION 

Geometrical frustration with strong electronic corre- 
lations is one of the main issues in modern condensed- 
matter physics. In connection with recent experimen- 
tal studies of frustrated quantum magnets such as those 
on triangular, kagome, spinel and pyrochlore lattices [1- 
7] as well as on triangular structure of 3 He on graphite 
[8], Hubbard model on lattices with geometrically frus- 
trated structures have been intensively studied [9-22]. 
For example, a quantum spin-liquid feature in the or- 
ganic material k-(ET)2Cu2(CN)3 stimulated the studies 
on the Hubbard model at half-filling on anisotropic trian- 
gular lattices [11, 13-15, 17, 21]. Another example is the 
two-dimensional Hubbard model with the next-nearest- 
neighbor hopping, t' , which will be studied in this paper. 
Although geometrical frustration tends to suppress con- 
ventional magnetic orderings, the properties of the resul- 
tant quantum phases and their competitions with other 
phases are, in many cases, still under discussion. 
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FIG. 1: The Lattice structure (a) and the hopping integrals 
t, t' and t± (a)(b) used in this study. 



In this paper, we study electronic properties of the 
Hubbard model at half-filling on the stacked square lat- 
tice as shown in Fig. 1. The Hamiltonian is given by 

H = -£( ]C * c ti,«r c i.».» - J2 ^ 4,l,<r c i^) 



-fi}n it i t(7 + U/nj,i,tTii,i,i- 



(1) 
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Here I labels the layers. We consider three different hop- 
pings t, t' and t± where t (t 1 ) is nearest- (next-nearest-) 
neighbor hopping in the same layer, and t± is interlayer 
hopping. Note that signs of t, t' and t± do not affect 
the results due to particle-hole symmetry for half-filled 
case. The operator c^i^ {c ilrJ ) annihilates (creates) an 
electron with spin a at the site i in the layer I, and riij,, a 
is the corresponding occupation number operator. \i is 
the chemical potential, and U is the local Coulomb re- 
pulsion. At half-filling, the supression of conventional 
magnetic orders and their competitions are still under 
discussion. 

In the large U region, the half-filled Hubbard model 
in Eq. (1) with t± — can be transformed into the 
two-dimensional J1-J2 Heisenberg model on the square 
lattice. This J\-Ji model has been studied extensively 
[23-35]. When J 2 = or t'/t = 0, it is well accepted that 
the ground state exhibits antiferromagnetic long-range 
order with the magnetic wave vector Qo = (71", it) due to 
a perfect nesting of fermi surface. As J2 increases, the an- 
tiferromagnetism is increasingly frustrated. In the limit 
of J2/J1 = 00, a collinear magnetic state with the mag- 
netic wave vector Qi = (it, 0) or (0,7r) is realized. This 
state consists of two antiferromagnetic long-range orders 
formed on the two sublattices. In the small U region, 
Mizusaki et al. recently studied the two-dimensional 



Hubbard model, Eq. (1) with t± = 0, using path-integral 
rcnormalization group method [21]. They found a long- 
period antiferromagnetic- insulator phase with 2x4 struc- 
ture for intermediate t'/t~0.7 at large U/t > 7. Fur- 
thermore, a quantum spin liquid phase with gapless spin 
excitations is found near the Mott transition to paramag- 
netic metals. They also claim that there is degeneracy of 
the ground states with various total momenta within the 
whole Brillouin zone in the quantum spin liquid phase. 
Therefore, various unconventional states are expected in 
the frustrated Hubbard model. 

The main problem we would like to study in this paper 
is the influence of frustration and dimensionality, namely 
that of t'/t and t±/t, on the magnetic property of the 
frustrated Hubbard model in Eq. (1). First, we will show 
that a paramagnetic insulating state, similar to that ob- 
tained by Mizusaki et al., is obtained in a different nu- 
merical method for 0.71 < t'/t < 0.85 and U/t = 10. 
Then we will study its stability by increasing the three 
dimensionality (ij_). Actually, the dimensionality also 
has strong influence on the property of magnetism [36] . 
Although the tendency to order is more pronounced in 
three-dimensional systems than in low-dimensional ones, 
a magnetically disordered phase can be seen in frustrated 
systems such as a pyrochlore and a stacked kagome lat- 
tice [37, 38]. Therefore, high-dimensionalization (adding 
t±/t), which naively seems to have the competing effect 
with the effect of t'/t on the ordering phenomenon, will 
affect the ground-state property. 

The exact diagonalization used in the two-dimensional 
Hubbard model will not be appropriate for the three- 
dimensional problem under consideration. Therefore, 
we use a quantum cluster approach [39] to describe 
strong electron correlations and geometrical frustration. 
In recent years there has been a substantial progress 
in numerical techniques using quantum cluster theories, 
such as cluster extensions of dynamical mean-field the- 
ory (DMFT) [40], i.e., dynamical cluster approximation 
(DC A) [41] and cellular DMFT (CDMFT) [42], or varia- 
tional cluster approach (VCA) [43, 44]. In this paper, we 
apply VCA which was proposed recently. It is based on a 
sclf-energy-functional theory (SFT) [45] which provides 
a general variational scheme using dynamical informa- 
tion from an exactly solvable "reference system" (in the 
present case, an isolated cluster) in order to study the 
infinite-size lattice fermion problem. It has been recently 
applied to the broken symmetry phases in the Hubbard 
model [46—49] and the dimer Hubbard model as a model 
for the organic materials k-(ET)2X [15]. 

This paper is organized as follows: We start with a 
brief review of SFT in general (Sec II A), and then de- 
scribe VCA applied in the Hamiltonian (1) (Sec II B). 
Some technical details are addressed in Appendix. In Sec. 
Ill, our results (energies, order parameters and a phase 
diagram) are presented and discussed. Finally, Sec. IV 
contains our main conclusions and a summary. 



II. THEORETICAL BACKGROUND 

In order to describe the VCA used in this paper, it is 
useful to explain the SFT first. 

A. Self-energy-functional theory[45] 

For a system with Hamiltonian H = Ho(t) + i?i(U), 
where t represents the one-particle and U the interaction 
parameters, the grand potential of the system at temper- 
ature, T, and chemical potential, /i, can be written as a 
functional of the self-energy S: 



fi t ,u[Z] = Tr ln(G^ - S)" 1 + Fu[E], 



(2) 



with the stational property <Jfit,u Pphys] = f° r the 
physical self-energy. Here, Go,t = (w + n — t) _1 
is the free Green's function and Tr is defined as 



Tr 



TJ2 U 



o iVn0 + 



tr, with tr being the usual trace and 



uj n being the Matsubara frequencies, u>„ = (2n+l)7rT, for 
integer n. Fjj[S] is the Legendre transform of the univer- 
sal Luttinger-Ward Functional $u[G] [52]. Because the 
Luttingcr-Ward Functional is defined via an infinite sum 
of renormalized skeleton diagrams, the functional depen- 
dence $u[G] or i*u[S] is not known explicitly. However, 
it is important that the functional form i<u[S] is inde- 
pendent of t . 

Due to this universality of Fu [S] , we have 

£V,u[S] - Tr ln(G-j, - S)" 1 + F V [H], (3) 

for the self-energy functional of a reference system, which 
is given by a Hamiltonian with the same interaction part 
U but modified one-particle parameters t', i.e. H' = 
Ho(t') + -ffi(U). By a proper choice of Ho(t'), the prob- 
lem posed by the reference system H ' can be much sim- 
pler than the original probrem posed by H , such that the 
self-energy of the reference system Et',u can be com- 
puted exactly within a certain subspace of parameters 
t' which we call a subspace, S. Combining Eqs. (2) and 
(3), we can eliminate the functional i*\j[E]. Inserting the 
self-energy of the reference system as a trial self-energy, 
we obtain 

fit,upV,u] = ^t',u + Tr ln(G -j - St^u)- 1 

-Tr lnGt'.u, (4) 

where f2t',u an d Gf,u = (GjT t , — S t '.u) _1 are the grand 
potential and the Green's function of the reference sys- 
tem. Stationaly points are obtained in the restricted sub- 
space S of trial self-energy (Et',u € S). Varying the 
trial self-energy in S is carried out by varying the one- 
particle parameters t' of the reference system. In addi- 
tion to the variational parameters, t', suitably chosen fic- 
titious symmetry-breaking Weiss fields can be introduced 
as variational parameters. By this method, normal and 
off-diagonal long-range order can be studied. Note that 
CPT, VCA and CDMFT, which treat decoupled clusters 
differently, are understood from a unified framework of 
SFT. 



B. Variational cluster approach 

VCA is conceptually clear and simple, and much easier 
to implement numerically. Therefore, we apply VCA to 
the Hamiltonian Eq. (1). We rewrite Eq. (1) as 



H = H (t)+H 1 , 

where Hi is the interaction part, 



H X = U 



id 



(5) 



(G) 



and -ffo(t) is the rest of the Hamiltonian Eq. (1). 

In the following, we fix the parameter, U, as U/t = 10, 
and the average particle number as (n^i^) + (^u,i) = lj 
and consider at zero temperature. The energy scale is set 
by choosing t = 1. 

We assume that the Hamiltonian of the reference sys- 
tem H' is given by a set of decoupled clusters of a finite 
size. For an individual cluster, the Hamiltonian reads 



H' 

cluster 



-^Hub 



#AF 



-^CM 



^local' 



(7) 



Here, H^ uh is the Hubbard Hamiltonian in the cluster, 
H' AF and H' CM are two symmetry-breaking terms (Weiss 
fields), 



#af 



/i'a„> Wt - ^,u)e jQAF ' Rl 



'AF, 



(8) 




FIG. 2: Superclusters which tile the three-dimensional lat- 
tice (a) with antiferromagnetic long-range order and (b) with 
collinear magnetic long-range order. Each supercluster con- 
sists of two 2x2 clusters which stack in the z- direct ion. A 
major spin at each site is shown by an blue arrow. 



Here 6(w) is the Heaviside step function and oj' m are 
the poles of Gt' : u, i-e., the one-particle excitation ener- 
gies of the cluster, u>' m = E r — E s , which is obtained by 
exact diagonalization. Here, we introduce the notation, 
to = (r, s), to indicate an excitaion between two states s 
and r. Similarly, uo m are the poles of the Green's function 
(Gpj - Et'.u) -1 - In E qs- (11) and (12), R represents 
a contribution due to the poles of the self-energy which 
cancels out in Eq. (4) and can thus be ignored. Numeri- 
cal techniques to calculate u) m and the Green's function 
of the supercluster are addressed in Appendixes A andB. 



and 



Hqm 



h'cM22( n i,l,1 



,i)e 



'iQcivi ■!*•;,, 



(9) 



where h' AF and h' CM are the strength of the fields used as 
variational parameters. Qaf = (tt, 7t, 7t) is the antiferro- 
magnetic (AF) wave vector, and Qcm = (vr, 0, n) is the 
collinear magnetic (CM) wave vector. Furthermore, in 
order to ensure thermodynamic consistency with respect 
to the average particle number, the site-independent en- 
ergy, e', is also treated as a variational parameter [47] 
which appears in the local term 

tfl'ocal = e'E 71 ^- ( 10 ) 

We use a three-dimensional cluster (called as superclus- 
ter) with 2x2x2 sites which consists of two 2x2 clusters 
stacked in the z-direction as shown in Fig. 2. and diag- 
onalize the Hamiltonian, H' cluater . 

The trace in Eq. (4) is evaluated accurately as follows. 
By converting the frequency integrals to a sum over the 
poles of the Green's function [45], we have 



Tr ln(Gol - Sf.u)" 



=°E^ 



,e(- 



and 



Tr InGv.u T =° J>™ (- 



R. 



R (11) 



(12) 



III. RESULTS 
A. Total energy 

First, let us consider the behavior of the total energy 
per site Eq. Figure 3 shows the obtained values of Eq 
for ij_ = as a function of t' . In this case, layers are 
indepenedent of each other. Eq's for paramagnetic Mott 
insulating (MI) state without magnetic long-range order, 
AF and CM states are denoted by solid, dashed and dash- 
dotted lines, respectively. Note that these three states 
are insulating because of gap opening in the total den- 
sity of states (not shown here). Figure 3 indicates that 
three different phases appear as a ground state for differ- 
ent values of t', namely that two magnetically long-range 
ordered states appear at small and large t' They are sep- 
arated by MI state in the parameter region t' c <t'< t' , 
where t^wO.71 and t^wO.85. E for MI (solid line) and 
Eq for AF (dashed line) smoothly coinside as increas- 
ing t' , indicating the second-order phase transition. On 
the other hand, E for MI (solid line) and Eq for CM 
(dash-dotted line) cross each other. Furthermore, hys- 
teresis behavior is observed at the transition point i C2 , 
indicating that the transition from MI to CM is of the 
first order. 

Figure 4 shows the same type of data for a rather 
strong interlayer coupling, t±_ = 0.7 which indicates 
that only two different phases appear as a ground state, 



namely AF state for t' < t' and CM state for t' > t' C3 , 
respectively, with t' c «0.77. In this case MI state does 
not appear and the transition from AF to CM is of the 
first order. The comparison of Figs. 3 and 4 shows that 
the MI state, which is stable in a rather wide t' region 
in the strictly two-dimensional system, becomes unstable 
by including the three-dimensionality. 
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t'/t 

FIG. 3: t' dependence of the total energy per site, Eo, for 
t± = obtained in three different states: paramagnetic Mott 
insulating (MI) (solid line), antiferromagnetic (AF) (dashed 
line) and collinear magnetic (CM) states (dash-dotted line). 
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FIG. 4: Same as in Fig. 3, for t± = 0.7. 



B. Order parameters and phase diagram 

Figure 5 shows the t' dependences of the AF order 
parameter at small t' and CM order parameter at large 
t' . In the region where no order parameter appears, MI 
state is stable as a ground state. As expected, the or- 
der parameters are monotonously increasing with t±_ and 
the phase transition points t' Cl , t' C2 changes as a func- 
tion of t±, which is shown in Fig. 6. At i^«0.44, the 



transition points t' ci and t' C2 coinside. For t± < tj_, AF 
order parameter vanishes continuously as is typical for 
the second-order transition. On the other hand, CM or- 
der parameter abruptly drops to zero, which is consistent 
with the first-order transition. For tj_ > t* ± , the values of 
AF order parameter and CM order parameter are discon- 
tinuously connected, indicating the first-order transition. 
The obtained phase diagram is summarised in Fig. 6, 
where solid line denotes the second-order transition and 
dashed line denotes the first-order transition. 




FIG. 5: t' dependence of antiferromagnetic (AF) and 
collinear magnetic (CM) order parameters for various val- 
ues of t±. At small t' , AF order parameters for t±/t = 
0, 0.3, 0.4, 0.5, 0.6, 0.7 are shown from bottom to top. At large 
t' , CM order parameters for the same values of t± are shown 
from bottom to top. 



Here let us compare our results with t±_ — with those 
in the J\-Ji Heisenberg model on the square lattice. It 
has been shown that for J2/J1 < 0.4, the AF state ap- 
pears and for J2/J1 > 0.6, the CM state appears [33]. 
Since J 1 ^^- and J2~-^j-, these critical values corre- 
spond to t'/i«0.63 and 0.77 which roughly agree with 
t' cl and t' c2 obtained in the present method. For the in- 
termediate region of 0.4 < J2/J1 < 0.6, no definite con- 
clusion has been drawn on the nature of the ground state 
of J1-J2 Heisenberg model. Possibilities of columnar- 
dimerized state [31, 53, 54], plaquette singlet state [55], 
and resonating-valence-bond state have been discussed. 
The intermediate state obtained in our method is a non- 
magnetic insulator. The identification of this state is a 
remaining future problem. 

Schmalfufi et.al. studied the stacked 3\-J% Heisen- 
berg model using a coupled-cluster method and rotation- 
invariant Green's function method [35]. They showed 
similar phase diagrams as ours, indicating the desappear- 
ance of an intermediate magnetically disordered quan- 
tum paramagnetic phase for quite small J±/Ji~0.2 — 0.3. 
This critical value corresponds to ij_/i«0.45 — 0.55 which 
agree again with Fig. 6. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
t'/t 

FIG. 6: Phase diagram in the t'-tx plane. Transition from 
antiferromagnetic (AF) state to paramagnetic Mott insulat- 
ing (MI) state is of the second-order (denoted by solid line). 
Transition from collinear magnetic (CM) state to MI or AF 
state is of the first-order (denoted by dashed line) . 



IV. SUMMARY AND CONCLUSIONS 



To summarize, we have studied the role of frustra- 
tion and dimensionality in the t-t' Hubbard model on 
the stacked square lattice given by the Hamiltonian, Eq. 
(1), at half-filling and at zero temparature, by changing 
t' and t±. We have employed recently proposed varia- 
tional cluster approach (VCA). Important advantages of 
VCA are that local and off-site short-range correlations 
are captured exactly, and that broken symmetries are 
treated with a rigorous dynamical variational principle. 
Therefore, VCA is a powerful method for analyzing mag- 
netic properties of Hubbard model with strong electron 
correlations and geometrical frustration. 

Our results show that t' destroys AF long-range order- 
ing continuously at small t' and induces CM long-range 
ordering at large t' with the first-order transition. There 
is a MI state at intermediate t' . We also find that t± 
has a tendency to stabilize the magnetic long-range or- 
ders and finally MI state vanishes at a critical value of t± 
(£^_«0.44). Therefore, interesting phenomena induced by 
strong electron correlations and geometrical frustration 
are very sensitive to dimensionality. In other words, ge- 
ometrical frustration and low dimensionality have a key 
role of intriguing physics at least in our model in this 
study. 
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APPENDIX A: CALCULATION OF THE POLES 
OF THE LATTICE GREEN'S FUNCTION 

The poles uj m of the Green's function (G^J — 
Xt^u) -1 = Gt,u can be obtained in the following way 
[48]: Consider the Lehmann representation of Gf,u 
which can be cast into the form 



Go,fl,t',u(w) — y Qa.m ~rQr, 

*■ — * i ,) — t, v 



m,/3' 



(Al) 



where a = (site i, spin a) in the cluster. The "Q-matrix" 
is defined as 



T=0 



= S r>0 (0\c a \s) + £ 8j0 (r|c a |0), 



(A2) 



and uj' m and Q a , m are obtained from exact diagonal- 
ization of the cluster. The spectral weight (residue) of 
Ga i/3 , t ',uM at a pole u> = u' m is given by Q a , m Q m ^- |0) 
denotes the (grand-canonical) ground state of the refer- 
ence system. Introducing the diagonal matrix g m . n (uj) — 

6 m ,n/(u ~ u'm)> we have: 



g*,u(w) = QsMQ f - 



(A3) 



Defining V — t — t', the VCA expression for the lattice 
Green's function can be written as 



Gt,u= 



1 



1 



G 0.t _ S t',U 



Gvu-V 



(A4) 



This expression can be transformed with the help of the 
Q-matrix in Eqs. (A2) and (A3): 

1 „ 1 



G t 



u 



(QgQ 1 ) 



Q- 



g~ 



QtVQ 



Q f . (A5) 



S m , n u' m , the poles of 



Since g _1 = uj — A with A mr 
Gt,u are now simply given by the eigenvalues of the (fre- 
quency independent) matrix M = A + Q T VQ and can 
be easily found by numerical diagonalization. The di- 
mension of M is given by the number of poles of Gt',u 
with nonvanishing spectral weight. 



APPENDIX B: CALCULATION OF THE 
SUPERCLUSTER GREEN'S FUNCTION 

The Green's function of the supercluster can be calcu- 
lated as follows: Switching off the hopping processes that 



connect the N c = 4 clusters in the supercluster gives a 
block-diagonal Hamiltonian which can be treated by nu- 
merical diagonalization. The switched off hopping pro- 
cesses are then incorporated again perturbatively. The 
supercluster Green's function G^/'u can be written as 

1 



W'U 



„ (block)- 1 
U t',U 



Y(s.c) 



(Bl) 



where 



^ (block) 
W'.U 



G 



(i) 



t',U 





-(2) 



QWg^MQWt 


Q« 
Q( 2 ) 

Q«t 
QWt 







g (1) M o 
o g ( 2 )(c) 



(2)t 



_ n (block) (block) 



Q 



g v 



>(w)Q 



(block)f 



(B2) 



Here, superscripts (1) and (2) denote a cluster index in a 
supercluster. V( sc ' is the intercluster hopping between 
the two clusters. The expression Eq. (Bl) can be trans- 
formed with the help of the Q( block )-cxpression of GL u° 
of Eq. (B2): 



i(s.c) _ 



G(S.C.) _ 
t',u - 



1 



CQ(block)„(block)Q(block)t')-l _ y(s.c) 
__ q (block) ' 



r (block)-l 
x Q (block) t_ 



Q(block)fy(s.c.)Q (block) 

(B3) 



Since g ( bl ° ck )-! = w - A( block ) where A( block ) is the 

diagonal matrix the elements of which are the poles 

of g^'{u>) and g^ 2 '(w), the poles of G^'u are given 

by the eigenvalues of the matrix M^ sc ' = A^ block ^ + 
Q(biock)t v (s. c .)Q(biock) and can be found by numerica i 

diagonalization. Diagonalizing MA S ' C ') by an appropriate 
matrix P as 



M (s.c.) = p M ( s - c )pt ) 



(B4) 



F (s.c) 



where M is the diagonal matrix the elements of which 
are the poles of Gif^j , Gi^ can be expressed as 



i(s.c) 



G;;;u J H = Q lsx - , r-'HQ 



,(s.c.)t 



(B5) 



with 



g 



(B.C.) 



M 



1 



Q(S.C.) 
Q( S .C.)t 



bj — M 

q (block) p 

p tQ(block)t^ 



(B6) 

(B7) 
(B8) 



Since Eq. (B5) has the same form as Eq. (A3), one 
repeat the procedure from Eq. (A3) to Eq. (A5) once 
again in order to obtain the poles of the lattice Green's 
function. 
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